Math
using System;
using System.Collections.Generic;
using System.Text;
namespace algorithms.math
{
// ----- Big (Long Arithmetics) --------------------------------------------
// Big()
// Big(string s)
// Big(int n), n < RADIX
// Big(Big big)
// int Length
// string ToString()
// void Abs()
// int CompareTo(Big big)
// Big Add(Big big)
// Big Sub(Big big)
// Big Mul(int k), k < RADIX
// Big Quotient(int k), k < RADIX
// int Remainder(int k), k < RADIX
// Big Power(int k), 0 <= k < RADIX
// Big Mul(Big big)
// static int Compare(Big big1, Big big2)
// static Big Add(Big big1, Big big2)
// static Big Sub(Big big1, Big big2)
// static Big Mul(Big big, int k), k < RADIX
// static Big Quotient(Big big, int k), k < RADIX
// static int Remainder(Big big, int k), k < RADIX
// static Big Power(Big big, int k), 0 <= k < RADIX
// static Big Mul(Big big1, Big big2)
// static Big Fact(int n), n < RADIX
// -------------------------------------------------------------------------
public class Big
{
const int RADIX = 1000 * 1000 * 1000;
const int RADIX_WIDTH = 9;
bool Neg { get; set; }
List<int> Num { get; set; }
public Big()
{
Neg = false;
Num = new List<int> { 0 };
}
public Big(string s)
{
Neg = s[0] == '-';
Num = new List<int>();
if (Neg) s = s.Substring(1);
s = s.TrimStart('0');
if (string.IsNullOrEmpty(s)) s = "0";
int end = s.Length;
while (end > 0)
{
int w = Math.Min(RADIX_WIDTH, end);
Num.Add(int.Parse(s.Substring(end - w, w)));
end -= w;
}
}
public Big(int n)
{
Neg = n < 0;
Num = new List<int>() { n };
}
public Big(Big big)
{
Neg = big.Neg;
Num = new List<int>(big.Num);
}
public int Length { get { return Num.Count; } }
public override string ToString()
{
string s = NumToString(Num);
if (Neg) s = "-" + s;
return s;
}
public void Abs()
{
Neg = false;
}
public int CompareTo(Big big)
{
return Compare(this, big);
}
public Big Add(Big big)
{
return Add(this, big);
}
public Big Sub(Big big)
{
return Sub(this, big);
}
public Big Mul(int k)
{
return Mul(this, k);
}
public Big Quotient(int k)
{
return Quotient(this, k);
}
public int Remainder(int k)
{
return Remainder(this, k);
}
public Big Power(int k)
{
return Power(this, k);
}
public Big Mul(Big big)
{
return Mul(this, big);
}
public static int Compare(Big big1, Big big2)
{
if (big1.Neg != big2.Neg) return big1.Neg ? -1 : 1;
if (big1.Neg) return Compare(big2.Num, big1.Num);
return Compare(big1.Num, big2.Num);
}
static Big AddSub(Big big1, Big big2, bool add)
{
if (add)
return new Big() { Neg = big1.Neg, Num = Big.Sum(big1.Num, big2.Num) };
else
{
int cmp = Big.Compare(big1.Num, big2.Num);
switch (cmp)
{
case -1:
return new Big() { Neg = !big1.Neg, Num = Big.Dif(big2.Num, big1.Num) };
case 1:
return new Big() { Neg = big1.Neg, Num = Big.Dif(big1.Num, big2.Num) };
}
return new Big(0);
}
}
public static Big Add(Big big1, Big big2)
{
return AddSub(big1, big2, big1.Neg == big2.Neg);
}
public static Big Sub(Big big1, Big big2)
{
return AddSub(big1, big2, big1.Neg != big2.Neg);
}
public static Big Mul(Big big, int k)
{
return new Big() { Neg = big.Neg != (k < 0), Num = Mul(big.Num, Math.Abs(k)) };
}
public static Big Quotient(Big big, int k)
{
return new Big() { Neg = big.Neg != (k < 0), Num = Quotient(big.Num, Math.Abs(k)) };
}
public static int Remainder(Big big, int k)
{
return Remainder(big.Num, Math.Abs(k));
}
public static Big Power(Big big, int k)
{
return new Big() { Neg = big.Neg && (k % 2 != 0), Num = Power(big.Num, k) };
}
public static Big Mul(Big big1, Big big2)
{
return new Big() { Neg = big1.Neg != big2.Neg, Num = Mul(big1.Num, big2.Num) };
}
public static Big Fact(int n)
{
Big big = new Big(1);
for (int i = 2; i <= n; i++) big = Mul(big, i);
return big;
}
static string NumToString(List<int> num)
{
StringBuilder sb = new StringBuilder();
foreach (int b in num)
sb.Insert(0, b.ToString().PadLeft(RADIX_WIDTH, '0'));
string s = sb.ToString().TrimStart('0');
if (string.IsNullOrEmpty(s)) s = "0";
return s;
}
static int Compare(List<int> num1, List<int> num2)
{
Norm(num1);
Norm(num2);
int cmp = num1.Count.CompareTo(num2.Count);
int ix = num1.Count;
while (cmp == 0 && ix > 0)
{
ix--;
cmp = num1[ix].CompareTo(num2[ix]);
}
return cmp;
}
static List<int> Sum(List<int> num1, List<int> num2)
{
int len = Math.Max(num1.Count, num2.Count);
List<int> sum = new List<int>(new int[len]);
int ix = 0;
bool carry = false;
while (ix < len || carry)
{
if (ix == sum.Count) sum.Add(0);
int b1 = ix < num1.Count ? num1[ix] : 0;
int b2 = ix < num2.Count ? num2[ix] : 0;
int b = b1 + b2 + (carry ? 1 : 0);
carry = b >= RADIX;
if (carry) b -= RADIX;
sum[ix] = b;
ix++;
}
return sum;
}
static List<int> Dif(List<int> num1, List<int> num2)
{
List<int> dif = new List<int>(new int[num1.Count]);
int ix = 0;
bool carry = false;
while (ix < num2.Count || carry)
{
int b1 = num1[ix];
int b2 = ix < num2.Count ? num2[ix] : 0;
int b = b1 - b2 - (carry ? 1 : 0);
carry = b < 0;
if (carry) b += RADIX;
dif[ix] = b;
ix++;
}
while (ix < num1.Count)
{
dif[ix] = num1[ix];
ix++;
}
Norm(dif);
return dif;
}
static List<int> Mul(List<int> num, int k)
{
List<int> product = new List<int>(new int[num.Count]);
int ix = 0;
int carry = 0;
while (ix < num.Count || carry > 0)
{
if (ix == product.Count) product.Add(0);
long nx = ix < num.Count ? num[ix] : 0;
long cur = nx * k + carry;
product[ix] = (int)(cur % RADIX);
carry = (int)(cur / RADIX);
ix++;
}
Norm(product);
return product;
}
static List<int> Quotient(List<int> num, int k)
{
List<int> quotient = new List<int>(new int[num.Count]);
int carry = 0;
for (int ix = num.Count - 1; ix >= 0; ix--)
{
long cur = num[ix] + (long)carry * RADIX;
quotient[ix] = (int)(cur / k);
carry = (int)(cur % k);
}
Norm(quotient);
return quotient;
}
static int Remainder(List<int> num, int k)
{
int carry = 0;
for (int ix = num.Count - 1; ix >= 0; ix--)
{
carry = (int)((num[ix] + (long)carry * RADIX) % k);
}
return carry;
}
static List<int> Power(List<int> num, int k)
{
List<int> power = new List<int>(new int[] { 1 });
while (k > 0)
{
if ((k & 1) == 1) power = Mul(power, num);
num = Mul(num, num);
k >>= 1;
}
Norm(power);
return power;
}
static List<int> Mul(List<int> num1, List<int> num2)
{
List<int> c = new List<int>(new int[num1.Count + num2.Count]);
for (int ix = 0; ix < num1.Count; ix++)
{
int iy = 0;
int carry = 0;
while (iy < num2.Count || carry > 0)
{
int b = iy < num2.Count ? num2[iy] : 0;
long cur = c[ix + iy] + (long)num1[ix] * b + carry;
c[ix + iy] = (int)(cur % RADIX);
carry = (int)(cur / RADIX);
iy++;
}
}
Norm(c);
return c;
}
static void Norm(List<int> num)
{
while (num.Count > 1 && num[num.Count - 1] == 0) num.RemoveAt(num.Count - 1);
}
}
// -------------------------------------------------------------------------
}
namespace algorithms.math
{
public static class BitHacks
{
// ----- Bit Hacks -----------------------------------------------------
//
// http://grepcode.com/file/repository.grepcode.com/java/root/jdk/openjdk/6-b14/java/lang/Long.java
//
// ulong highestOneBit(ulong i)
// long lowestOneBit(long i)
// int numberOfLeadingZeros(long i)
// int numberOfTrailingZeros(long i)
// int bitCount(long i)
// long rotateLeft(long i, int distance)
// long rotateRight(long i, int distance)
// long reverse(long i)
// long reverseBytes(long i)
//
// bool is_a_power_of_2(uint v)
// uint lexicographically_next_bit_permutation(uint v)
// ---------------------------------------------------------------------
public static ulong highestOneBit(ulong i)
{
i |= (i >> 1);
i |= (i >> 2);
i |= (i >> 4);
i |= (i >> 8);
i |= (i >> 16);
i |= (i >> 32);
return i - (i >> 1);
}
public static long lowestOneBit(long i)
{
return i & -i;
}
public static int numberOfLeadingZeros(long i)
{
if (i == 0) return 64;
int n = 1;
uint x = (uint)((ulong)i >> 32);
if (x == 0) { n += 32; x = (uint)i; }
if (x >> 16 == 0) { n += 16; x <<= 16; }
if (x >> 24 == 0) { n += 8; x <<= 8; }
if (x >> 28 == 0) { n += 4; x <<= 4; }
if (x >> 30 == 0) { n += 2; x <<= 2; }
n -= (int)(x >> 31);
return n;
}
public static int numberOfTrailingZeros(long i)
{
uint x, y;
if (i == 0) return 64;
int n = 63;
y = (uint)i; if (y != 0) { n = n - 32; x = y; } else x = (uint)(i >> 32);
y = x << 16; if (y != 0) { n = n - 16; x = y; }
y = x << 8; if (y != 0) { n = n - 8; x = y; }
y = x << 4; if (y != 0) { n = n - 4; x = y; }
y = x << 2; if (y != 0) { n = n - 2; x = y; }
return n - (int)((x << 1) >> 31);
}
public static int bitCount(long i)
{
i = i - (long)(((ulong)i >> 1) & 0x5555555555555555L);
i = (i & 0x3333333333333333L) + (long)(((ulong)i >> 2) & 0x3333333333333333L);
i = (i + (long)((ulong)i >> 4)) & 0x0f0f0f0f0f0f0f0fL;
i = i + (long)((ulong)i >> 8);
i = i + (long)((ulong)i >> 16);
i = i + (long)((ulong)i >> 32);
return (int)i & 0x7f;
}
public static long rotateLeft(long i, int distance)
{
return (i << distance) | (long)((ulong)i >> -distance);
}
public static long rotateRight(long i, int distance)
{
return (long)((ulong)i >> distance) | (i << -distance);
}
public static long reverse(long i)
{
i = (i & 0x5555555555555555L) << 1 | (long)((ulong)i >> 1) & 0x5555555555555555L;
i = (i & 0x3333333333333333L) << 2 | (long)((ulong)i >> 2) & 0x3333333333333333L;
i = (i & 0x0f0f0f0f0f0f0f0fL) << 4 | (long)((ulong)i >> 4) & 0x0f0f0f0f0f0f0f0fL;
i = (i & 0x00ff00ff00ff00ffL) << 8 | (long)((ulong)i >> 8) & 0x00ff00ff00ff00ffL;
i = (i << 48) | ((i & 0xffff0000L) << 16) | (long)(((ulong)i >> 16) & 0xffff0000L) | (long)((ulong)i >> 48);
return i;
}
public static long reverseBytes(long i)
{
i = (i & 0x00ff00ff00ff00ffL) << 8 | (long)((ulong)i >> 8) & 0x00ff00ff00ff00ffL;
return (i << 48) | ((i & 0xffff0000L) << 16) | (long)(((ulong)i >> 16) & 0xffff0000L) | (long)((ulong)i >> 48);
}
public static bool is_a_power_of_2(uint v)
{
return (v & (v - 1)) == 0;
}
public static uint lexicographically_next_bit_permutation(uint v)
{
uint t = (v | (v - 1)) + 1;
return t | (uint)((((t & -t) / (v & -v)) >> 1) - 1);
}
// ---------------------------------------------------------------------
}
}
namespace algorithms.math
{
public static class Bits
{
// ----- Bits ----------------------------------------------------------
//
// int[] BitsArray(int nbits)
// void MarkBit(int[] bits, int bit)
// void ClearBit(int[] bits, int bit)
// bool IsMarked(int[] bits, int bit)
//
// void MarkBit(ref int bits, int bit)
// void ClearBit(ref int bits, int bit)
// bool IsMarked(ref int bits, int bit)
// ---------------------------------------------------------------------
public static int[] BitsArray(int nbits)
{
return new int[(nbits - 1) / 32 + 1];
}
public static void MarkBit(int[] bits, int bit)
{
bits[bit / 32] |= (1 << (bit % 32));
}
public static void ClearBit(int[] bits, int bit)
{
bits[bit / 32] &= ~(1 << (bit % 32));
}
public static bool IsMarked(int[] bits, int bit)
{
return (bits[bit / 32] & (1 << (bit % 32))) != 0;
}
public static void MarkBit(ref int bits, int bit)
{
bits |= (1 << bit);
}
public static void ClearBit(ref int bits, int bit)
{
bits &= ~(1 << bit);
}
public static bool IsMarked(ref int bits, int bit)
{
return (bits & (1 << bit)) != 0;
}
// ---------------------------------------------------------------------
}
}
namespace algorithms.math
{
public static class Bits32
{
// ----- Bits32 --------------------------------------------------------
//
// int MarkBit(int bits, int bit)
// int ClearBit(int bits, int bit)
// bool IsMarked(int bits, int bit)
// uint RemoveBit(uint bits, int bit)
// ---------------------------------------------------------------------
public static int MarkBit(int bits, int bit)
{
return bits | (1 << bit);
}
public static int ClearBit(int bits, int bit)
{
return bits & ~(1 << bit);
}
public static bool IsMarked(int bits, int bit)
{
return (bits & (1 << bit)) != 0;
}
static uint[] B32 = new uint[] {
0x1, 0x2, 0x4, 0x8,
0x10, 0x20, 0x40, 0x80,
0x100, 0x200, 0x400, 0x800,
0x1000, 0x2000, 0x4000, 0x8000,
0x10000, 0x20000, 0x40000, 0x80000,
0x100000, 0x200000, 0x400000, 0x800000,
0x1000000, 0x2000000, 0x4000000, 0x8000000,
0x10000000, 0x20000000, 0x40000000, 0x80000000,
};
static uint[] F32 = new uint[] {
0x0,
0x1, 0x3, 0x7, 0xF,
0x1F, 0x3F, 0x7F, 0xFF,
0x1FF, 0x3FF, 0x7FF, 0xFFF,
0x1FFF, 0x3FFF, 0x7FFF, 0xFFFF,
0x1FFFF, 0x3FFFF, 0x7FFFF, 0xFFFFF,
0x1FFFFF, 0x3FFFFF, 0x7FFFFF, 0xFFFFFF,
0x1FFFFFF, 0x3FFFFFF, 0x7FFFFFF, 0xFFFFFFF,
0x1FFFFFFF, 0x3FFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF
};
static uint RemoveBit(uint bits, int bit)
{
return ((bits >> 1) & (~F32[bit])) | (bits & F32[bit]);
}
static uint ReverseBits(uint bits, int n)
{
uint revs = 0;
for (int i = 0; i < n; i++)
if ((bits & B32[i]) == B32[i]) revs |= B32[n - 1 - i];
return revs;
}
// ---------------------------------------------------------------------
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace algorithms.math
{
public static class Combinatorics
{
// ----- Combinatorics -------------------------------------------------
//
// IEnumerable<T[]> AllFor<T>(T[] array)
// IEnumerable<int[]> Combinations(int m, int n)
// https://en.wikipedia.org/wiki/Heap%27s_algorithm
// void Heap<T>(T[] A, Action doPermutation)
// ---------------------------------------------------------------------
public static IEnumerable<T[]> Permutations<T>(T[] array)
{
if (array == null || array.Length == 0)
{
yield return new T[0];
}
else
{
for (int pick = 0; pick < array.Length; ++pick)
{
T item = array[pick];
int i = -1;
T[] rest = System.Array.FindAll<T>(
array, delegate (T p) { return ++i != pick; }
);
foreach (T[] restPermuted in Permutations(rest))
{
i = -1;
yield return System.Array.ConvertAll<T, T>(
array,
delegate (T p)
{
return ++i == 0 ? item : restPermuted[i - 1];
}
);
}
}
}
}
public static IEnumerable<int[]> Combinations(int m, int n)
{
int[] result = new int[m];
Stack<int> stack = new Stack<int>();
stack.Push(0);
while (stack.Count > 0)
{
int index = stack.Count - 1;
int value = stack.Pop();
while (value < n)
{
result[index++] = ++value;
stack.Push(value);
if (index == m)
{
yield return result;
break;
}
}
}
}
public static void Heap<T>(T[] A, Action doPermutation)
{
int N = A.Length;
int[] C = new int[N];
doPermutation();
int i = 0;
while (i < N)
if (C[i] < i)
{
int j = i % 2 == 0 ? 0 : C[i];
T t = A[j]; A[j] = A[i]; A[i] = t;
doPermutation();
C[i]++;
i = 0;
}
else
{
C[i] = 0;
i++;
}
}
// ---------------------------------------------------------------------
}
}
namespace algorithms.math
{
public static class Common
{
// ----- Common --------------------------------------------------------
//
// int gcd(int a, int b)
// ---------------------------------------------------------------------
public static int gcd(int a, int b)
{
while (b != 0) b = a % (a = b);
return a;
}
// ---------------------------------------------------------------------
}
}
using System;
namespace algorithms.math
{
// ----- Complex -----------------------------------------------------------
// double Re
// double Im
// Complex(double re, double im)
// Complex(double re)
// double Abs()
// double Abs2()
// Complex Copy()
// double Arg()
// Complex Conjugate()
// Complex Unit()
// int CompareTo(object o)
// bool Equals(object o)
// int GetHashCode()
// string ToString()
// static Complex Zero
// static Complex I
// static Complex MaxValue
// static Complex MinValue
// static Complex FromAbsArg(double abs, double arg)
// static explicit operator Complex(double f)
// static explicit operator double(Complex c)
// static bool operator ==(Complex a, Complex b)
// static bool operator !=(Complex a, Complex b)
// static Complex operator +(Complex a)
// static Complex operator -(Complex a)
// static Complex operator +(Complex a, double f)
// static Complex operator +(double f, Complex a)
// static Complex operator +(Complex a, Complex b)
// static Complex operator -(Complex a, double f)
// static Complex operator -(double f, Complex a)
// static Complex operator -(Complex a, Complex b)
// static Complex operator *(Complex a, double f)
// static Complex operator *(double f, Complex a)
// static Complex operator *(Complex a, Complex b)
// static Complex operator /(Complex a, double f)
// static Complex operator /(Complex a, Complex b)
// static bool IsEqual(Complex a, Complex b, double tolerance)
// -------------------------------------------------------------------------
public struct Complex: IComparable
{
public double Re { get; private set; }
public double Im { get; private set; }
public Complex(double re, double im) { Re = re; Im = im; }
public Complex(double re) { Re = re; Im = 0; }
public double Abs() { return Math.Sqrt(Re * Re + Im * Im); }
public double Abs2() { return Re * Re + Im * Im; }
public Complex Copy() { return new Complex(Re, Im); }
public double Arg() { return Math.Atan2(Im, Re); }
public Complex Conjugate() { return new Complex(Re, -Im); }
public Complex Unit() { return this / Abs(); }
public int CompareTo(object o) { if (o is Complex) { return Abs().CompareTo(((Complex)o).Abs()); } if (o is double) { return Abs().CompareTo((double)o); } return 0; }
public override bool Equals(object o) { if (o is Complex) { Complex c = (Complex)o; return (this == c); } return false; }
public override int GetHashCode() { return (Re.GetHashCode() ^ Im.GetHashCode()); }
public override string ToString() { return string.Format("({0}, {1})", Re, Im); }
public static Complex Zero { get { return new Complex(0, 0); } }
public static Complex I { get { return new Complex(0, 1); } }
public static Complex MaxValue { get { return new Complex(double.MaxValue, double.MaxValue); } }
public static Complex MinValue { get { return new Complex(double.MinValue, double.MinValue); } }
public static Complex FromAbsArg(double abs, double arg) { return new Complex(abs * Math.Cos(arg), abs * Math.Sin(arg)); }
public static explicit operator Complex(double f) { return new Complex(f, 0); }
public static explicit operator double(Complex c) { return c.Re; }
public static bool operator ==(Complex a, Complex b) { return (a.Re == b.Re) && (a.Im == b.Im); }
public static bool operator !=(Complex a, Complex b) { return (a.Re != b.Re) || (a.Im != b.Im); }
public static Complex operator +(Complex a) { return new Complex(a.Re, a.Im); }
public static Complex operator -(Complex a) { return new Complex(-a.Re, -a.Im); }
public static Complex operator +(Complex a, double f) { return new Complex(a.Re + f, a.Im); }
public static Complex operator +(double f, Complex a) { return new Complex(a.Re + f, a.Im); }
public static Complex operator +(Complex a, Complex b) { return new Complex(a.Re + b.Re, a.Im + b.Im); }
public static Complex operator -(Complex a, double f) { return new Complex(a.Re - f, a.Im); }
public static Complex operator -(double f, Complex a) { return new Complex(f - a.Re, a.Im); }
public static Complex operator -(Complex a, Complex b) { return new Complex(a.Re - b.Re, a.Im - b.Im); }
public static Complex operator *(Complex a, double f) { return new Complex(a.Re * f, a.Im * f); }
public static Complex operator *(double f, Complex a) { return new Complex(a.Re * f, a.Im * f); }
public static Complex operator *(Complex a, Complex b) { return new Complex(a.Re * b.Re - a.Im * b.Im, a.Re * b.Im + a.Im * b.Re); }
public static Complex operator /(Complex a, double f) { return new Complex(a.Re / f, a.Im / f); }
public static Complex operator /(Complex a, Complex b) {
double denom = b.Re * b.Re + b.Im * b.Im;
return new Complex((a.Re * b.Re + a.Im * b.Im) / denom, (a.Im * b.Re - a.Re * b.Im) / denom);
}
public static bool IsEqual(Complex a, Complex b, double tolerance)
{
return (Math.Abs(a.Re - b.Re) < tolerance) && (Math.Abs(a.Im - b.Im) < tolerance);
}
}
// -------------------------------------------------------------------------
}
using System;
using System.Linq;
namespace algorithms.math
{
public static class FastFourierTransform
{
// ----- Fast Fourier Transform ----------------------------------------
//
// Uses Cooley-Tukey iterative in-place algorithm with radix-2 DIT case
// assumes no of points provided are a power of 2
//
// Depends on:
// -- Complex (algorithms.math)
//
// void FFT(Complex[] buffer, int n, bool invert = false)
// void FFT(Complex[] buffer, bool invert)
// int[] FFTMultiply(int[] a, int na, int[] b, int nb)
// int[] FFTMultiply(int[] a, int[] b)
// ---------------------------------------------------------------------
public static void FFT(Complex[] buffer, int n, bool invert = false)
{
int dig = 0;
while ((1 << dig) < n) dig++;
int[] rev = new int[n];
for (int i = 0; i < n; i++)
{
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (dig - 1));
if (rev[i] > i)
{
var temp = buffer[i];
buffer[i] = buffer[rev[i]];
buffer[rev[i]] = temp;
}
}
for (int len = 2; len <= n; len <<= 1)
{
double angle = 2 * Math.PI / len;
if (invert) angle *= -1;
Complex wgo = new Complex(Math.Cos(angle), Math.Sin(angle));
for (int i = 0; i < n; i += len)
{
Complex w = new Complex(1);
for (int j = 0; j < (len >> 1); j++)
{
Complex a = buffer[i + j];
Complex b = w * buffer[i + j + (len >> 1)];
buffer[i + j] = a + b;
buffer[i + j + (len >> 1)] = a - b;
w *= wgo;
}
}
}
if (invert)
for (int i = 0; i < n; i++) buffer[i] /= n;
}
public static void FFT(Complex[] buffer, bool invert = false)
{
FFT(buffer, buffer.Length, invert);
}
public static int[] FFTMultiply(int[] a, int na, int[] b, int nb)
{
int sz = Math.Max(na, nb);
int n = 1;
while (n < sz) n <<= 1;
n <<= 1;
Complex[] fa = new Complex[n];
for (int i = 0; i < na; i++) fa[i] = new Complex(a[i]);
for (int i = na; i < n; i++) fa[i] = Complex.Zero;
Complex[] fb = new Complex[n];
for (int i = 0; i < nb; i++) fb[i] = new Complex(b[i]);
for (int i = nb; i < n; i++) fb[i] = Complex.Zero;
FFT(fa);
FFT(fb);
for (int i = 0; i < n; i++) fa[i] *= fb[i];
FFT(fa, true);
return fa.Select(p => (int)(p.Re < 0 ? p.Re - 0.5 : p.Re + 0.5)).ToArray();
}
public static int[] FFTMultiply(int[] a, int[] b)
{
return FFTMultiply(a, a.Length, b, b.Length);
}
// ---------------------------------------------------------------------
}
}
using System;
using System.Numerics;
namespace algorithms.math
{
// ----- Frac (Fraction) ---------------------------------------------------
// Frac()
// Frac(long num)
// Frac(long num, long den)
// Frac(BigInteger num)
// Frac(BigInteger num, BigInteger den)
// Frac Reduce()
// static Frac Add(Frac a, Frac b)
// static Frac Sub(Frac a, Frac b)
// static Frac Mul(Frac a, Frac b)
// static Frac Div(Frac a, Frac b)
// static Frac Inv(Frac a)
// Frac Add(Frac b)
// Frac Sub(Frac b)
// Frac Mul(Frac b)
// Frac Div(Frac b)
// Frac Inv()
// string ToString()
// string ToStringSimple()
// int CompareTo(Frac f)
// bool SimpleEquals(Frac f)
// bool Equals(Object o)
// int GetHashCode()
// -------------------------------------------------------------------------
public class Frac: IComparable<Frac>
{
public BigInteger num { get; private set; }
public BigInteger den { get; private set; }
public Frac() { this.num = BigInteger.Zero; this.den = BigInteger.One; }
public Frac(long num) { this.num = new BigInteger(num); this.den = BigInteger.One; }
public Frac(long num, long den) { this.num = new BigInteger(num); this.den = new BigInteger(den); Reduce(); }
public Frac(BigInteger num) { this.num = num; this.den = BigInteger.One; }
public Frac(BigInteger num, BigInteger den) { this.num = num; this.den = den; Reduce(); }
public Frac Reduce()
{
if (den.Sign == 0) { }
else if (num.Sign == 0) { den = BigInteger.One; }
else
{
if (den.Sign < 0)
{
num = BigInteger.Negate(num);
den = BigInteger.Negate(den);
}
BigInteger g = BigInteger.GreatestCommonDivisor(num, den);
num = BigInteger.Divide(num, g);
den = BigInteger.Divide(den, g);
}
return this;
}
public static Frac Add(Frac a, Frac b) { return new Frac(BigInteger.Add(BigInteger.Multiply(a.num, b.den), BigInteger.Multiply(a.den, b.num)), BigInteger.Multiply(a.den, b.den)); }
public static Frac Sub(Frac a, Frac b) { return new Frac(BigInteger.Subtract(BigInteger.Multiply(a.num, b.den), BigInteger.Multiply(a.den, b.num)), BigInteger.Multiply(a.den, b.den)); }
public static Frac Mul(Frac a, Frac b) { return new Frac(BigInteger.Multiply(a.num, b.num), BigInteger.Multiply(a.den, b.den)); }
public static Frac Div(Frac a, Frac b) { return new Frac(BigInteger.Multiply(a.num, b.den), BigInteger.Multiply(a.den, b.num)); }
public static Frac Inv(Frac a) { return new Frac(a.den, a.num); }
public Frac Add(Frac b) { num = BigInteger.Add(BigInteger.Multiply(num, b.den), BigInteger.Multiply(den, b.num)); den = BigInteger.Multiply(den, b.den); return Reduce(); }
public Frac Sub(Frac b) { num = BigInteger.Subtract(BigInteger.Multiply(num, b.den), BigInteger.Multiply(den, b.num)); den = BigInteger.Multiply(den, b.den); return Reduce(); }
public Frac Mul(Frac b) { num = BigInteger.Multiply(num, b.num); den = BigInteger.Multiply(den, b.den); return Reduce(); }
public Frac Div(Frac b) { num = BigInteger.Multiply(num, b.den); den = BigInteger.Multiply(den, b.num); return Reduce(); }
public Frac Inv() { BigInteger d = num; num = den; den = d; return Reduce(); }
public override string ToString() { return num.ToString() + "/" + den.ToString(); }
public string ToStringSimple() { return den.Equals(BigInteger.One) ? num.ToString() : num.ToString() + "/" + den.ToString(); }
public int CompareTo(Frac f) { return BigInteger.Multiply(num, f.den).CompareTo(BigInteger.Multiply(f.num, den)); }
public bool SimpleEquals(Frac f) { return num.Equals(f.num) && den.Equals(f.den); }
public override bool Equals(Object o)
{
if (o == null) return false;
Frac f = (Frac)o;
return BigInteger.Multiply(num, f.den).Equals(BigInteger.Multiply(f.num, den));
}
public override int GetHashCode() { return (num.GetHashCode() ^ den.GetHashCode()); }
}
// -------------------------------------------------------------------------
}
using System;
namespace algorithms.math
{
public static class KthElementInArray
{
// ----- Kth Element In Array ------------------------------------------
//
// A randomized algorithm based on QuickSort partition
// ASSUMPTION: ELEMENTS IN A[] ARE DISTINCT
//
// O(n)
//
// int RandomizedSelect<T>(T[] a, int k)
// int RandomizedSelect<T>(T[] a, int ix, int len, int k)
// ---------------------------------------------------------------------
static Random rnd = new Random();
public static int RandomizedSelect<T>(T[] a, int k) where T: IComparable
{
return RandomizedSelect(a, 0, a.Length, k);
}
public static int RandomizedSelect<T>(T[] a, int ix, int len, int k) where T : IComparable
{
if (len == 1) return ix;
int q = RandomizedPartition(a, ix, len);
int m = q - ix;
if (k < m) return RandomizedSelect(a, ix, m, k);
return RandomizedSelect(a, q, len - m, k - m);
}
static int RandomizedPartition<T>(T[] a, int ix, int len) where T : IComparable
{
int i = rnd.Next(len);
T temp = a[ix];
a[ix] = a[ix + i];
a[ix + i] = temp;
return Partition(a, ix, len);
}
static int Partition<T>(T[]a, int ix, int len) where T : IComparable
{
T x = a[ix];
int i = ix;
int j = ix + len - 1;
while (true)
{
while (a[j].CompareTo(x) >= 0) { j--; }
while (a[i].CompareTo(x) < 0) { i++; }
if (i >= j) return i;
T temp = a[i];
a[i] = a[j];
a[j] = temp;
}
}
// ---------------------------------------------------------------------
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace algorithms.math
{
// ----- Matrix ------------------------------------------------------------
// int Rows
// int Cols
// Matrix(int rows, int cols)
// bool IsSquare()
// double this[int row, int col]
// Matrix GetCol(int k)
// void SetCol(Matrix v, int k)
// void MakeLU()
// Matrix SolveWith(Matrix v)
// void MakeRref()
// Matrix Invert()
// double Det()
// Matrix GetP()
// Matrix Duplicate()
// string ToString()
// static Matrix SubsForth(Matrix A, Matrix b)
// static Matrix SubsBack(Matrix A, Matrix b)
// static Matrix ZeroMatrix(int rows, int cols)
// static Matrix IdentityMatrix(int rows, int cols)
// static Matrix RandomMatrix(int rows, int cols, int dispersion)
// static Matrix Parse(string ps)
// static Matrix Transpose(Matrix m)
// static Matrix Power(Matrix m, int pow)
// static Matrix operator -(Matrix m)
// static Matrix operator +(Matrix m1, Matrix m2)
// static Matrix operator -(Matrix m1, Matrix m2)
// static Matrix operator *(Matrix m1, Matrix m2)
// static Matrix operator *(double n, Matrix m)
// -------------------------------------------------------------------------
public class Matrix
{
public int Rows { get; private set; }
public int Cols { get; private set; }
private double[] mat;
private Matrix L;
private Matrix U;
private int[] pi;
private double detOfP = 1;
public Matrix(int rows, int cols)
{
Rows = rows;
Cols = cols;
mat = new double[rows * cols];
}
public bool IsSquare()
{
return (Rows == Cols);
}
public double this[int row, int col]
{
get { return mat[row * Cols + col]; }
set { mat[row * Cols + col] = value; }
}
public Matrix GetCol(int k)
{
Matrix m = new Matrix(Rows, 1);
for (int i = 0; i < Rows; i++) m[i, 0] = this[i, k];
return m;
}
public void SetCol(Matrix v, int k)
{
for (int i = 0; i < Rows; i++) this[i, k] = v[i, 0];
}
public void MakeLU()
{
if (!IsSquare()) throw new MException("The matrix is not square!");
L = IdentityMatrix(Rows, Cols);
U = Duplicate();
pi = new int[Rows];
for (int i = 0; i < Rows; i++) pi[i] = i;
double p = 0;
double pom2;
int k0 = 0;
int pom1 = 0;
for (int k = 0; k < Cols - 1; k++)
{
p = 0;
// find the row with the biggest pivot
for (int i = k; i < Rows; i++)
{
if (Math.Abs(U[i, k]) > p)
{
p = Math.Abs(U[i, k]);
k0 = i;
}
}
if (p == 0) throw new MException("The matrix is singular!");
// switch two rows in permutation matrix
pom1 = pi[k]; pi[k] = pi[k0]; pi[k0] = pom1;
for (int i = 0; i < k; i++)
{
pom2 = L[k, i]; L[k, i] = L[k0, i]; L[k0, i] = pom2;
}
if (k != k0) detOfP *= -1;
// Switch rows in U
for (int i = 0; i < Cols; i++)
{
pom2 = U[k, i]; U[k, i] = U[k0, i]; U[k0, i] = pom2;
}
for (int i = k + 1; i < Rows; i++)
{
L[i, k] = U[i, k] / U[k, k];
for (int j = k; j < Cols; j++)
U[i, j] = U[i, j] - L[i, k] * U[k, j];
}
}
}
// Function solves Ax = v in confirmity with solution vector "v"
public Matrix SolveWith(Matrix v)
{
if (Rows != Cols) throw new MException("The matrix is not square!");
if (Rows != v.Rows) throw new MException("Wrong number of results in solution vector!");
if (L == null) MakeLU();
Matrix b = new Matrix(Rows, 1);
// switch two items in "v" due to permutation matrix
for (int i = 0; i < Rows; i++) b[i, 0] = v[pi[i], 0];
Matrix z = SubsForth(L, b);
Matrix x = SubsBack(U, z);
return x;
}
// Function makes reduced echolon form
public void MakeRref()
{
int lead = 0;
for (int r = 0; r < Rows; r++)
{
if (Cols <= lead) break;
int i = r;
while (this[i, lead] == 0)
{
i++;
if (i == Rows)
{
i = r;
lead++;
if (Cols == lead)
{
lead--;
break;
}
}
}
for (int j = 0; j < Cols; j++)
{
double temp = this[r, j];
this[r, j] = this[i, j];
this[i, j] = temp;
}
double div = this[r, lead];
for (int j = 0; j < Cols; j++) this[r, j] /= div;
for (int j = 0; j < Rows; j++)
{
if (j != r)
{
double sub = this[j, lead];
for (int k = 0; k < Cols; k++) this[j, k] -= (sub * this[r, k]);
}
}
lead++;
}
}
// Function returns the inverted matrix
public Matrix Invert()
{
if (L == null) MakeLU();
Matrix inv = new Matrix(Rows, Cols);
for (int i = 0; i < Rows; i++)
{
Matrix Ei = Matrix.ZeroMatrix(Rows, 1);
Ei[i, 0] = 1;
Matrix col = SolveWith(Ei);
inv.SetCol(col, i);
}
return inv;
}
// Function for determinant
public double Det()
{
if (L == null) MakeLU();
double det = detOfP;
for (int i = 0; i < Rows; i++) det *= U[i, i];
return det;
}
// Function returns permutation matrix "P" due to permutation vector "pi"
public Matrix GetP()
{
if (L == null) MakeLU();
Matrix matrix = ZeroMatrix(Rows, Cols);
for (int i = 0; i < Rows; i++) matrix[pi[i], i] = 1;
return matrix;
}
// Function returns the copy of this matrix
public Matrix Duplicate()
{
Matrix matrix = new Matrix(Rows, Cols);
for (int i = 0; i < Rows; i++)
for (int j = 0; j < Cols; j++)
matrix[i, j] = this[i, j];
return matrix;
}
// Function solves Ax = b for A as a lower triangular matrix
public static Matrix SubsForth(Matrix A, Matrix b)
{
if (A.L == null) A.MakeLU();
int n = A.Rows;
Matrix x = new Matrix(n, 1);
for (int i = 0; i < n; i++)
{
x[i, 0] = b[i, 0];
for (int j = 0; j < i; j++) x[i, 0] -= A[i, j] * x[j, 0];
x[i, 0] = x[i, 0] / A[i, i];
}
return x;
}
// Function solves Ax = b for A as an upper triangular matrix
public static Matrix SubsBack(Matrix A, Matrix b)
{
if (A.L == null) A.MakeLU();
int n = A.Rows;
Matrix x = new Matrix(n, 1);
for (int i = n - 1; i > -1; i--)
{
x[i, 0] = b[i, 0];
for (int j = n - 1; j > i; j--) x[i, 0] -= A[i, j] * x[j, 0];
x[i, 0] = x[i, 0] / A[i, i];
}
return x;
}
// Function generates the zero matrix
public static Matrix ZeroMatrix(int rows, int cols)
{
Matrix matrix = new Matrix(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
matrix[i, j] = 0;
return matrix;
}
// Function generates the identity matrix
public static Matrix IdentityMatrix(int rows, int cols)
{
Matrix matrix = ZeroMatrix(rows, cols);
for (int i = 0; i < Math.Min(rows, cols); i++)
matrix[i, i] = 1;
return matrix;
}
// Function generates the random matrix
public static Matrix RandomMatrix(int rows, int cols, int dispersion)
{
Random random = new Random();
Matrix matrix = new Matrix(rows, cols);
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
matrix[i, j] = random.Next(-dispersion, dispersion);
return matrix;
}
// Function parses the matrix from string
public static Matrix Parse(string ps)
{
string s = NormalizeMatrixString(ps);
string[] rows = System.Text.RegularExpressions.Regex.Split(s, "\r\n");
string[] nums = rows[0].Split(' ');
Matrix matrix = new Matrix(rows.Length, nums.Length);
try
{
for (int i = 0; i < rows.Length; i++)
{
nums = rows[i].Split(' ');
for (int j = 0; j < nums.Length; j++) matrix[i, j] = double.Parse(nums[j]);
}
}
catch (FormatException) { throw new MException("Wrong input format!"); }
return matrix;
}
// Function returns matrix as a string
public override string ToString()
{
StringBuilder s = new StringBuilder();
for (int i = 0; i < Rows; i++)
{
for (int j = 0; j < Cols; j++)
s.Append(String.Format("{0,5:E2}", this[i, j]) + " ");
s.AppendLine();
}
return s.ToString();
}
// Matrix transpose, for any rectangular matrix
public static Matrix Transpose(Matrix m)
{
Matrix t = new Matrix(m.Cols, m.Rows);
for (int i = 0; i < m.Rows; i++)
for (int j = 0; j < m.Cols; j++)
t[j, i] = m[i, j];
return t;
}
// Power matrix to exponent
public static Matrix Power(Matrix m, int pow)
{
if (pow == 0) return IdentityMatrix(m.Rows, m.Cols);
if (pow == 1) return m.Duplicate();
if (pow == -1) return m.Invert();
Matrix x;
if (pow < 0) { x = m.Invert(); pow *= -1; }
else x = m.Duplicate();
Matrix ret = IdentityMatrix(m.Rows, m.Cols);
while (pow != 0)
{
if ((pow & 1) == 1) ret *= x;
x *= x;
pow >>= 1;
}
return ret;
}
private static void SafeAplusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
{
for (int i = 0; i < size; i++) // rows
for (int j = 0; j < size; j++) // cols
{
C[i, j] = 0;
if (xa + j < A.Cols && ya + i < A.Rows) C[i, j] += A[ya + i, xa + j];
if (xb + j < B.Cols && yb + i < B.Rows) C[i, j] += B[yb + i, xb + j];
}
}
private static void SafeAminusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
{
for (int i = 0; i < size; i++) // rows
for (int j = 0; j < size; j++) // cols
{
C[i, j] = 0;
if (xa + j < A.Cols && ya + i < A.Rows) C[i, j] += A[ya + i, xa + j];
if (xb + j < B.Cols && yb + i < B.Rows) C[i, j] -= B[yb + i, xb + j];
}
}
private static void SafeACopytoC(Matrix A, int xa, int ya, Matrix C, int size)
{
for (int i = 0; i < size; i++) // rows
for (int j = 0; j < size; j++) // cols
{
C[i, j] = 0;
if (xa + j < A.Cols && ya + i < A.Rows) C[i, j] += A[ya + i, xa + j];
}
}
private static void AplusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
{
for (int i = 0; i < size; i++) // rows
for (int j = 0; j < size; j++) C[i, j] = A[ya + i, xa + j] + B[yb + i, xb + j];
}
private static void AminusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
{
for (int i = 0; i < size; i++) // rows
for (int j = 0; j < size; j++) C[i, j] = A[ya + i, xa + j] - B[yb + i, xb + j];
}
private static void ACopytoC(Matrix A, int xa, int ya, Matrix C, int size)
{
for (int i = 0; i < size; i++) // rows
for (int j = 0; j < size; j++) C[i, j] = A[ya + i, xa + j];
}
// Smart matrix multiplication
private static Matrix StrassenMultiply(Matrix A, Matrix B)
{
if (A.Cols != B.Rows) throw new MException("Wrong dimension of matrix!");
Matrix R;
int msize = Math.Max(Math.Max(A.Rows, A.Cols), Math.Max(B.Rows, B.Cols));
int size = 1; int n = 0;
while (msize > size) { size *= 2; n++; };
int h = size / 2;
Matrix[,] mField = new Matrix[n, 9];
/*
* 8x8, 8x8, 8x8, ...
* 4x4, 4x4, 4x4, ...
* 2x2, 2x2, 2x2, ...
* . . .
*/
int z;
for (int i = 0; i < n - 4; i++) // rows
{
z = (int)Math.Pow(2, n - i - 1);
for (int j = 0; j < 9; j++) mField[i, j] = new Matrix(z, z);
}
SafeAplusBintoC(A, 0, 0, A, h, h, mField[0, 0], h);
SafeAplusBintoC(B, 0, 0, B, h, h, mField[0, 1], h);
StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 1], 1, mField); // (A11 + A22) * (B11 + B22);
SafeAplusBintoC(A, 0, h, A, h, h, mField[0, 0], h);
SafeACopytoC(B, 0, 0, mField[0, 1], h);
StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 2], 1, mField); // (A21 + A22) * B11;
SafeACopytoC(A, 0, 0, mField[0, 0], h);
SafeAminusBintoC(B, h, 0, B, h, h, mField[0, 1], h);
StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 3], 1, mField); //A11 * (B12 - B22);
SafeACopytoC(A, h, h, mField[0, 0], h);
SafeAminusBintoC(B, 0, h, B, 0, 0, mField[0, 1], h);
StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 4], 1, mField); //A22 * (B21 - B11);
SafeAplusBintoC(A, 0, 0, A, h, 0, mField[0, 0], h);
SafeACopytoC(B, h, h, mField[0, 1], h);
StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 5], 1, mField); //(A11 + A12) * B22;
SafeAminusBintoC(A, 0, h, A, 0, 0, mField[0, 0], h);
SafeAplusBintoC(B, 0, 0, B, h, 0, mField[0, 1], h);
StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 6], 1, mField); //(A21 - A11) * (B11 + B12);
SafeAminusBintoC(A, h, 0, A, h, h, mField[0, 0], h);
SafeAplusBintoC(B, 0, h, B, h, h, mField[0, 1], h);
StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 7], 1, mField); // (A12 - A22) * (B21 + B22);
R = new Matrix(A.Rows, B.Cols); // result
/// C11
for (int i = 0; i < Math.Min(h, R.Rows); i++) // rows
for (int j = 0; j < Math.Min(h, R.Cols); j++) // cols
R[i, j] = mField[0, 1 + 1][i, j] + mField[0, 1 + 4][i, j] - mField[0, 1 + 5][i, j] + mField[0, 1 + 7][i, j];
/// C12
for (int i = 0; i < Math.Min(h, R.Rows); i++) // rows
for (int j = h; j < Math.Min(2 * h, R.Cols); j++) // cols
R[i, j] = mField[0, 1 + 3][i, j - h] + mField[0, 1 + 5][i, j - h];
/// C21
for (int i = h; i < Math.Min(2 * h, R.Rows); i++) // rows
for (int j = 0; j < Math.Min(h, R.Cols); j++) // cols
R[i, j] = mField[0, 1 + 2][i - h, j] + mField[0, 1 + 4][i - h, j];
/// C22
for (int i = h; i < Math.Min(2 * h, R.Rows); i++) // rows
for (int j = h; j < Math.Min(2 * h, R.Cols); j++) // cols
R[i, j] = mField[0, 1 + 1][i - h, j - h] - mField[0, 1 + 2][i - h, j - h] + mField[0, 1 + 3][i - h, j - h] + mField[0, 1 + 6][i - h, j - h];
return R;
}
// A * B into C, level of recursion, matrix field
private static void StrassenMultiplyRun(Matrix A, Matrix B, Matrix C, int l, Matrix[,] f)
{
int size = A.Rows;
int h = size / 2;
AplusBintoC(A, 0, 0, A, h, h, f[l, 0], h);
AplusBintoC(B, 0, 0, B, h, h, f[l, 1], h);
StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 1], l + 1, f); // (A11 + A22) * (B11 + B22);
AplusBintoC(A, 0, h, A, h, h, f[l, 0], h);
ACopytoC(B, 0, 0, f[l, 1], h);
StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 2], l + 1, f); // (A21 + A22) * B11;
ACopytoC(A, 0, 0, f[l, 0], h);
AminusBintoC(B, h, 0, B, h, h, f[l, 1], h);
StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 3], l + 1, f); //A11 * (B12 - B22);
ACopytoC(A, h, h, f[l, 0], h);
AminusBintoC(B, 0, h, B, 0, 0, f[l, 1], h);
StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 4], l + 1, f); //A22 * (B21 - B11);
AplusBintoC(A, 0, 0, A, h, 0, f[l, 0], h);
ACopytoC(B, h, h, f[l, 1], h);
StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 5], l + 1, f); //(A11 + A12) * B22;
AminusBintoC(A, 0, h, A, 0, 0, f[l, 0], h);
AplusBintoC(B, 0, 0, B, h, 0, f[l, 1], h);
StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 6], l + 1, f); //(A21 - A11) * (B11 + B12);
AminusBintoC(A, h, 0, A, h, h, f[l, 0], h);
AplusBintoC(B, 0, h, B, h, h, f[l, 1], h);
StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 7], l + 1, f); // (A12 - A22) * (B21 + B22);
/// C11
for (int i = 0; i < h; i++) // rows
for (int j = 0; j < h; j++) // cols
C[i, j] = f[l, 1 + 1][i, j] + f[l, 1 + 4][i, j] - f[l, 1 + 5][i, j] + f[l, 1 + 7][i, j];
/// C12
for (int i = 0; i < h; i++) // rows
for (int j = h; j < size; j++) // cols
C[i, j] = f[l, 1 + 3][i, j - h] + f[l, 1 + 5][i, j - h];
/// C21
for (int i = h; i < size; i++) // rows
for (int j = 0; j < h; j++) // cols
C[i, j] = f[l, 1 + 2][i - h, j] + f[l, 1 + 4][i - h, j];
/// C22
for (int i = h; i < size; i++) // rows
for (int j = h; j < size; j++) // cols
C[i, j] = f[l, 1 + 1][i - h, j - h] - f[l, 1 + 2][i - h, j - h] + f[l, 1 + 3][i - h, j - h] + f[l, 1 + 6][i - h, j - h];
}
// Stupid matrix multiplication
private static Matrix StupidMultiply(Matrix m1, Matrix m2)
{
if (m1.Cols != m2.Rows) throw new MException("Wrong dimensions of matrix!");
Matrix result = ZeroMatrix(m1.Rows, m2.Cols);
for (int i = 0; i < result.Rows; i++)
for (int j = 0; j < result.Cols; j++)
for (int k = 0; k < m1.Cols; k++)
result[i, j] += m1[i, k] * m2[k, j];
return result;
}
// Matrix multiplication
private static Matrix Multiply(Matrix m1, Matrix m2)
{
if (m1.Cols != m2.Rows) throw new MException("Wrong dimension of matrix!");
int msize = Math.Max(Math.Max(m1.Rows, m1.Cols), Math.Max(m2.Rows, m2.Cols));
// stupid multiplication faster for small matrices
if (msize < 32)
{
return StupidMultiply(m1, m2);
}
// stupid multiplication faster for non square matrices
if (!m1.IsSquare() || !m2.IsSquare())
{
return StupidMultiply(m1, m2);
}
// Strassen multiplication is faster for large square matrix 2^N x 2^N
// NOTE because of previous checks msize == m1.cols == m1.rows == m2.cols == m2.cols
double exponent = Math.Log(msize) / Math.Log(2);
if (Math.Pow(2, exponent) == msize)
{
return StrassenMultiply(m1, m2);
}
else
{
return StupidMultiply(m1, m2);
}
}
// Multiplication by constant n
private static Matrix Multiply(double n, Matrix m)
{
Matrix r = new Matrix(m.Rows, m.Cols);
for (int i = 0; i < m.Rows; i++)
for (int j = 0; j < m.Cols; j++)
r[i, j] = m[i, j] * n;
return r;
}
private static Matrix Add(Matrix m1, Matrix m2)
{
if (m1.Rows != m2.Rows || m1.Cols != m2.Cols) throw new MException("Matrices must have the same dimensions!");
Matrix r = new Matrix(m1.Rows, m1.Cols);
for (int i = 0; i < r.Rows; i++)
for (int j = 0; j < r.Cols; j++)
r[i, j] = m1[i, j] + m2[i, j];
return r;
}
public static string NormalizeMatrixString(string matStr)
{
// Remove any multiple spaces
while (matStr.IndexOf(" ") != -1)
matStr = matStr.Replace(" ", " ");
// Remove any spaces before or after newlines
matStr = matStr.Replace(" \r\n", "\r\n");
matStr = matStr.Replace("\r\n ", "\r\n");
// If the data ends in a newline, remove the trailing newline.
// Make it easier by first replacing \r\n’s with |’s then
// restore the |’s with \r\n’s
matStr = matStr.Replace("\r\n", "|");
while (matStr.LastIndexOf("|") == (matStr.Length - 1))
matStr = matStr.Substring(0, matStr.Length - 1);
matStr = matStr.Replace("|", "\r\n");
return matStr.Trim();
}
public static Matrix operator -(Matrix m)
{ return Matrix.Multiply(-1, m); }
public static Matrix operator +(Matrix m1, Matrix m2)
{ return Matrix.Add(m1, m2); }
public static Matrix operator -(Matrix m1, Matrix m2)
{ return Matrix.Add(m1, -m2); }
public static Matrix operator *(Matrix m1, Matrix m2)
{ return Matrix.Multiply(m1, m2); }
public static Matrix operator *(double n, Matrix m)
{ return Matrix.Multiply(n, m); }
}
public class MException : Exception
{
public MException(string Message) : base(Message) { }
}
// -------------------------------------------------------------------------
}
namespace algorithms.math
{
public static class ModularArithmetics
{
// ----- Modular Arithmetics -------------------------------------------
//
// long FastFib(long n, int modulus)
// long Multiply(long a, long b, int modulus)
// long Divide(long a, long b, int modulus)
// long Power(long a, long b, int modulus)
// long Choose(long n, int r, int modulus)
// ---------------------------------------------------------------------
// Fast doubling algorithm
public static long FastFib(long n, int modulus)
{
int k = n > int.MaxValue ? 63 : 31;
long a = 0;
long b = 1;
for (int i = k; i >= 0; i--)
{
long d = a * (b * 2 - a + modulus);
long e = a * a + b * b;
a = d % modulus;
b = e % modulus;
if ((((ulong)n >> i) & 1) != 0)
{
long c = a + b;
a = b;
b = c % modulus;
}
}
return a;
}
public static long Multiply(long a, long b, int modulus)
{
return (a * b) % modulus;
}
public static long Divide(long a, long b, int modulus)
{
return (a * Inverse(b, modulus)) % modulus;
}
public static long Power(long a, long b, int modulus)
{
long res = 1;
while (b > 0)
{
if ((b & 1) == 1) res = (res * a) % modulus;
a = (a * a) % modulus;
b >>= 1;
}
return res;
}
public static long Choose(long n, int r, int modulus)
{
long choose = 1;
while (true)
{
if (r == 0) break;
long N = n % modulus;
int R = r % modulus;
if (N < R) return 0;
for (int i = 0; i < R; i++)
choose = (choose * (N - i)) % modulus;
long factR = 1;
for (int i = 0; i < R; i++)
factR = (factR * (i + 1)) % modulus;
choose = Divide(choose, factR, modulus);
if (choose < 0) choose += modulus;
n /= modulus;
r /= modulus;
}
return choose;
}
static long Inverse(long a, int modulus)
{
long b = modulus;
long p = 1, q = 0;
while (b > 0)
{
long c = a / b;
long d = a;
a = b;
b = d % b;
d = p;
p = q;
q = d - c * q;
}
return p < 0 ? p + modulus : p;
}
// ---------------------------------------------------------------------
}
}
namespace algorithms.math
{
public static class ModularCombinatorics
{
// ----- Modular Combinatorics -----------------------------------------
//
// int[][] FiF(int nmax, int modulus)
// long Choose(int n, int r, int[][] fif, int modulus)
// ---------------------------------------------------------------------
public static int[][] FiF(int nmax, int modulus)
{
int[][] fif = new int[2][];
fif[0] = new int[nmax + 1];
fif[0][0] = 1;
for (int i = 1; i <= nmax; i++) fif[0][i] = (int)(((long)fif[0][i - 1] * i) % modulus);
long a = fif[0][nmax];
long b = modulus;
long p = 1;
long q = 0;
while (b > 0)
{
long c = a / b;
long d = a;
a = b;
b = d % b;
d = p;
p = q;
q = d - c * q;
}
fif[1] = new int[nmax + 1];
fif[1][nmax] = (int)(p < 0 ? p + modulus : p);
for (int i = nmax - 1; i >= 0; i--) fif[1][i] = (int)(((long)fif[1][i + 1] * (i + 1)) % modulus);
return fif;
}
public static long Choose(int n, int r, int[][] fif, int modulus)
{
if (n < 0 || r < 0 || r > n) return 0;
long factn = fif[0][n];
long invFactr = fif[1][r];
long invFactnr = fif[1][n - r];
long c1 = ((long)fif[0][n] * fif[1][r]) % modulus;
long c2 = (c1 * fif[1][n - r]) % modulus;
return ((((long)fif[0][n] * fif[1][r]) % modulus) * fif[1][n - r]) % modulus;
}
// ---------------------------------------------------------------------
}
}
using System;
using System.Collections.Generic;
namespace algorithms.math
{
public static class NumberTheory
{
// ----- Number Theory -------------------------------------------------
//
// Depends on:
// -- ModularArithmetics (algorithms.math)
//
// List<int> Primes(int limit)
// Dictionary<int, int> PrimeFactors(long n, List<int> primes)
// long Phi(long n, List<int> primes)
// bool IsPrimitiveRoot(long prime, long root, List<int> primes)
// ---------------------------------------------------------------------
public static List<int> Primes(int limit)
{
List<int> primes = new List<int>();
for (int k = 2; k <= limit; k++)
{
bool prime = true;
int rk = (int)Math.Sqrt(k);
foreach (int p in primes)
{
if (rk < p) break;
if (k % p == 0)
{
prime = false;
break;
}
}
if (prime) primes.Add(k);
}
return primes;
}
public static Dictionary<int, int> PrimeFactors(long n, List<int> primes)
{
Dictionary<int, int> pk = new Dictionary<int, int>();
int nr = (int)Math.Sqrt(n);
foreach (int p in primes)
{
if (n == 1 || p > nr) break;
if (n % p == 0)
{
pk[p] = 0;
while (n % p == 0)
{
pk[p]++;
n /= p;
}
}
}
if (n > 1) pk[(int)n] = 1;
return pk;
}
public static long Phi(long n, List<int> primes)
{
long phi = n;
foreach (int p in PrimeFactors(n, primes).Keys)
phi = phi / p * (p - 1);
return phi;
}
public static bool IsPrimitiveRoot(long prime, long root, List<int> primes)
{
foreach (int p in PrimeFactors(prime - 1, primes).Keys)
if (ModularArithmetics.Power(root, prime / p, (int)prime) == 1) return false;
return true;
}
// ---------------------------------------------------------------------
static bool Suspect(long b, int t, long u, long n)
{
long prod = 1;
while (u > 0)
{
if ((u & 1) == 1) prod = (prod * b) % n;
b = (b * b) % n;
u /= 2;
}
if (prod == 1) return true;
for (int i = 0; i < t; i++)
{
if (prod == n - 1) return true;
prod = (prod * prod) % n;
}
return false;
}
// up to 2^32-1
public static bool IsPrime(long n)
{
long k = n - 1;
int t = 0;
while (k % 2 == 0) { t++; k /= 2; }
if (n > 2 && n % 2 == 0) return false;
if (n > 3 && n % 3 == 0) return false;
if (n > 5 && n % 5 == 0) return false;
if (n > 7 && n % 7 == 0) return false;
return Suspect(61, t, k, n) && Suspect(7, t, k, n) && Suspect(2, t, k, n);
}
// up to 3*10^14
public static bool IsPrime314(long n)
{
if (n < 20 && (n == 2 || n == 3 || n == 5 || n == 7 || n == 11 || n == 13 || n == 17 || n == 19)) return true;
long k = n - 1;
int t = 0;
while (k % 2 == 0) { t++; k /= 2; }
if (n > 2 && n % 2 == 0) return false;
if (n > 3 && n % 3 == 0) return false;
if (n > 5 && n % 5 == 0) return false;
if (n > 7 && n % 7 == 0) return false;
return Suspect(19, t, k, n) && Suspect(17, t, k, n) && Suspect(13, t, k, n) && Suspect(11, t, k, n) && Suspect(7, t, k, n) && Suspect(5, t, k, n) && Suspect(3, t, k, n) && Suspect(2, t, k, n);
}
// ---------------------------------------------------------------------
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace algorithms.math
{
// ----- Polynomial --------------------------------------------------------
//
// polynomial = sum c_i * x^e_i
//
// static Polynomial Build(params Term[] t)
// static Polynomial Build(params long[] c)
// long Substitute(int v)
// static Polynomial Add(params Polynomial[] a)
// static Polynomial Mul(Polynomial a, long v)
// Polynomial Mul(long v)
// static Polynomial Mul(Polynomial a, Polynomial b)
// static Polynomial Pow(Polynomial a, int n)
// string ToString()
// -------------------------------------------------------------------------
public class Polynomial
{
public struct Term
{
public long C { get; set; }
public int E { get; set; }
public Term(long c, int e) { C = c; E = e; }
}
public Term[] terms { get; private set; }
private Polynomial() : base() { }
public static Polynomial Build(params Term[] t)
{
Polynomial p = new Polynomial();
p.terms = t;
return p.Simplify();
}
// x^0, x^1, ...
public static Polynomial Build(params long[] c)
{
Polynomial p = new Polynomial();
int n = c.Length;
p.terms = new Term[n];
for (int i = 0; i < n; i++)
{
p.terms[i] = new Term(c[i], i);
}
return p.Simplify();
}
public long Substitute(int v)
{
long ret = 0;
foreach (Term t in terms)
{
long lret = t.C;
for (int j = 0; j < t.E; j++) lret *= v;
ret += lret;
}
return ret;
}
public static Polynomial Add(params Polynomial[] a)
{
Polynomial p = new Polynomial();
int len = 0;
foreach (Polynomial spm in a) len += spm.terms.Length;
p.terms = new Term[len];
int q = 0;
foreach (Polynomial spm in a)
foreach (Term t in spm.terms)
p.terms[q++] = new Term(t.C, t.E);
return p.Simplify();
}
public static Polynomial Mul(Polynomial a, long v)
{
Polynomial p = new Polynomial();
int m = a.terms.Length;
p.terms = new Term[m];
for (int i = 0; i < m; i++)
p.terms[i] = new Term(a.terms[i].C * v, a.terms[i].E);
return p.Simplify();
}
public Polynomial Mul(long v)
{
for (int i = 0; i < terms.Length; i++) terms[i].C *= v;
return v == 0 ? Simplify() : this;
}
public static Polynomial Mul(Polynomial a, Polynomial b)
{
Polynomial p = new Polynomial();
int m = a.terms.Length, n = b.terms.Length;
p.terms = new Term[m * n];
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
Term t = new Term(
a.terms[i].C * b.terms[j].C,
a.terms[i].E + b.terms[j].E
);
p.terms[i * n + j] = t;
}
}
return p.Simplify();
}
static int numberOfLeadingZeros(long i)
{
if (i == 0) return 64;
int n = 1;
uint x = (uint)((ulong)i >> 32);
if (x == 0) { n += 32; x = (uint)i; }
if (x >> 16 == 0) { n += 16; x <<= 16; }
if (x >> 24 == 0) { n += 8; x <<= 8; }
if (x >> 28 == 0) { n += 4; x <<= 4; }
if (x >> 30 == 0) { n += 2; x <<= 2; }
n -= (int)(x >> 31);
return n;
}
public static Polynomial Pow(Polynomial a, int n)
{
if (a.terms.Length == 0) return a;
Polynomial ret = new Polynomial();
ret.terms = new Term[] { new Term(1L, 0) };
int x = 63 - numberOfLeadingZeros(n);
for (; x >= 0; x--)
{
ret = Mul(ret, ret);
if (n << 63 - x < 0) ret = Mul(ret, a);
}
return ret;
}
public override string ToString()
{
StringBuilder sb = new StringBuilder();
bool first = true;
foreach (Term t in terms)
{
if (t.C >= 0 && !first) sb.Append('+');
first = false;
if (t.C == -1 || t.C == 1)
{
if (t.E == 0)
sb.Append(t.C);
else
if (t.C == -1)
sb.Append('-');
}
else
sb.Append(t.C);
if (t.E > 0)
{
sb.Append('x');
if (t.E >= 2) sb.Append("^" + t.E);
}
}
return sb.ToString();
}
Polynomial Simplify()
{
if (terms.Length == 0) return this;
Array.Sort(terms, (a, b) => -(a.E - b.E));
int n = terms.Length;
int p = 1;
for (int i = 1; i < n; i++)
{
if (terms[i].E == terms[p - 1].E)
terms[p - 1].C += terms[i].C;
else
{
if (terms[p - 1].C == 0L) p--;
terms[p++] = terms[i];
}
}
if (terms[p - 1].C == 0L) p--;
if (terms.Length != p)
{
Term[] t = new Term[p];
Array.Copy(terms, t, p);
terms = t;
}
return this;
}
}
// -------------------------------------------------------------------------
}